% @author:夜剑听雨
% E-mail:2530595378@qq.com
% 2025.08.29
%%
% 频散数值模拟程序
clc;clear;close all;
%% 速度模型
%% 速度递增4层模型
% Vs = [150 300 400 600];           % 横波速度(m/s)
% Thickness = [2 4 6];              % 地层厚度(m)
%% 低速夹层4层模型
% Vs = [250 120 350 500];           % 横波速度(m/s)
% Thickness = [4 2 4];              % 地层厚度(m)

%% 高速夹层4层模型
Vs = [100 200 500 400];           % 横波速度(m/s)
Thickness = [4 2 4];              % 地层厚度(m)

%% 计算数值模型的频散曲线解析解，理论频散曲线
Vp = [1000 1000 1000 1000];       % 纵波速度(m/s)
Den = [2000 2000 2000 2000];      % 地层密度(kg/m3)
Freq = 5:65;                     % 频率范围
% 所有输入值必须是行向量
phase_velocity = calcmulti(Freq,Vs,Thickness,Vp,Den);
pv = phase_velocity(:,1);
%% 绘制理论频散曲线
figure;

% 非零元素索引
Ind = find(pv~=0);
% 绘制曲线
plot(Freq(Ind), pv(Ind), 'k', 'LineWidth', 1, ...
    'Marker', 'o', ...
    'MarkerSize', 3, ...
    'MarkerFaceColor', 'r');

% 设置y轴的范围和刻度
ylim([50 800]);
yticks(50:100:800);
% 设置x轴的范围和刻度
xlim([0 70]);
xticks(0:10:70);
% 设置图形的标题和轴标签
% 设置坐标轴线宽
set(gca, 'LineWidth', 1);
% Theoretical dispersion curves  'Model1
title('Theoretical dispersion curves','FontName', 'Microsoft YaHei', 'FontSize', 12);
xlabel('Frequency(Hz)','FontName', 'Microsoft YaHei', 'FontSize', 12);
ylabel('Phase velocity(m/s)','FontName', 'Microsoft YaHei', 'FontSize', 12);

% 消除图片白边
set(gca,'LooseInset',get(gca,'TightInset')+0.01);

% 关闭hold状态
hold off;

%% 台阶状速度图
% 画图最大深度
maxDepth = 20;
[vs_sta, depth_sta] = calcStairsData(Vs,Thickness, maxDepth);

% 创建图形
figure;
% 使用stairs函数绘制台阶状图形
stairs(vs_sta,depth_sta,  'r-', 'LineWidth', 1.5);
% % 设置y轴的范围和刻度
ylim([0 20]);
yticks(0:5:20);
% 设置x轴的范围和刻度
xlim([0 800]);
xticks(0:100:800);

% 设置x轴和y轴的标签
xlabel('S-wave velocity (m/s)','FontName', 'Microsoft YaHei', 'FontSize', 12);
ylabel('Depth (m)','FontName', 'Microsoft YaHei', 'FontSize', 12);
% 设置坐标轴线宽
set(gca, 'LineWidth', 1);
% 反转y轴
set(gca,'XAxisLocation','bottom');% 将x轴的位置设置在底部。
set(gca,'YDir','reverse'); % 将y轴方向设置为反向(从上到下递增)。
% 设置图形的标题速度递增模型
title('Velocity model','FontName', 'Microsoft YaHei', 'FontSize', 12);
% 消除图片白边
set(gca,'LooseInset',get(gca,'TightInset')+0.01);